function dif_B = PMCirFunc(r,theta,z,xd,yd,zd,flag)
%integrand of flux of a cylindrical PM
x = r.*cos(theta);
y = r.*sin(theta);
R1 = sqrt((xd-x).^2+(yd-y).^2+(zd-z).^2);

switch (flag)
    case 1 
        dif_B = r.*(xd-x)./R1.^3;
    case 2
        dif_B = r.*(yd-y)./R1.^3;
    case 3
        dif_B = r.*(zd-z)./R1.^3;
end
        